Structure Tensors

The first step to performing tractography research is in computing structure tensors. Structure tensors are essentially gradients displayed in matrix form. As a reminder, gradients themselves are essentially partial derivatives. We're interested in the gradient information because it's very useful for detecting certainty/uncertainty moving in directions in the neighborhood from a given point.

As discussed here, https://www.cs.cmu.edu/~sarsen/structureTensorTutorial/#7, if we have some structure tensor $S_w$:

As defined by the CAPTURE paper:

In the upper image, the p and r refer to voxel (spatial) locations. In the lower image, the $S_0$ matrix is the symmetric second-moment matrix generated from all the individual gradients ($I_x$ is the derivative w.r.t. x, for instance). By eigendecomposition of the $S_w$ matrix, we should get three eigenvectors (these can be combined to create a diffusion ellipsoid, which (if I'm understandng correctly) shows the certainty of movement in the three cardinal directions?).

Generating Structure Tensors

Following the MATLAB code provided by CAPTURE here: http://capture-clarity.org/clarity-based-tractography/

Most notably, in order to calculate the gradients, the MATLAB script uses an approach where it convolves the original image data with the derivative of 3D Gaussian kernels.

It then computed gradient products by taking the three gradients (in X, Y, Z, I'm guessing) and multiplied by themselves to get $I_x I_x$, $I_x I_y$, $I_x I_z$, $I_y I_y$, $I_y I_z$, $I_z I_z$. By the fact that the MATLAB script didn't compute $I_y I_x$, $I_z I_x$, and $I_z I_y$ I'm assuming they must be equal to the previously calculated values?

It then found the "amplitude" $g_a$ by finding: $\sqrt{I_x^2 + I_y^2 + I_z^2}$

It then saved this "amplitude" .nii. From there, it also saved a copy of the gradient vector $g_v$.

Now, to compute the structure tensor, the MATLAB CAPTURE script followed steps very similar to that of a difference of Gaussians technique. First the script generated a 3D Gaussian kernel with a given sigma value we defined at the beginning of the pipeline. From there, it convolved $I_x I_x$, $I_x I_y$, $I_x I_z$, $I_y I_y$, $I_y I_z$, $I_z I_z$ with the Gaussian kernel to obtain "blurred" versions of the original gradients. Then, these results were saved in two formats:

Structure Tensor Formats [in order of appearance in MATLAB concatenation]:
1) FSL format: blurred_XX, blurred_XY, blurred_XZ, blurred_YY, blurred_YZ, blurred_ZZ
2) DTK format: ? (was an odd combination of FSL, odd non-intuitive indexing)

These were saved as .nii files. Opening these files on SimpleITK showed:

Analysis:

1) This CAPTURE pipeline accepts TIF stacks, which I modified to accepting TIFF stacks. I used the 5 TIFF images from Aut1367 (the ones I used for Ilastik) as an example - these summed to 125 MB each, and the structure tensor .nii files from this were absolutely huge (3 GB for the FSL format, 850 MB for the DTK format).
2) Computationally, running these files on my computer for even such a small subset eventually froze my MATLAB. I have a hunch it's because MATLAB is storing all intermediary variables in memory. We'll need some way to deal with this either in MATLAB or in Python if we re-implement this.
3) The results: a lot of black. Intensities ranged from 0 to around 1.2 (surprisingly, they were floats... I was expecting integers). This may be because the sample data had no discernable tracts (so ask Greg on Thursday for suitable tractography data, and smaller data sets).
4) The next steps of the CAPTURE pipeline used the .nii files generated here and ran them through DiffusionToolkit to generate .trk files for visualization in TrackVis. However, the next steps require a "_brainmask.nii" and "_seedmask.nii" file to properly run the next steps. I'm unsure of what these are, and they haven't provided any clear examples.